Conway’s Game of Life

Game of life (code in order of presentation)

What are we building?

We are implementing GOL, resulting in something like the following

using Plots, Distributions, DynamicGrids

@time out = sim!(ArrayOutput(rand(Bernoulli(0.3), (100,100)), tspan=1:100), Life());
anim = @animate for i in 1:100
   heatmap(out[i], c=:binary, aspectratio=1) 
end;
gif(anim, fps=10)

  2.568484 seconds (6.82 M allocations: 374.577 MiB, 4.54% gc time, 99.55% compilation time)

┌ Info: Saved animation to 
│   fn = /home/michael/phd/GOL/tmp.gif
└ @ Plots /home/michael/.julia/packages/Plots/1KWPG/src/animation.jl:114

Starting by defining the rules of GOL

Rules

Any living cell with < 2 neighbors dies.
Any living cell with two or three live neighours lives
Any cell with > 3 neighbors dies
Any dead cell with 3 neighbors becomes alive

We want to design it not as a one-off simulation, by like we are designing a software package for GOL.

Why do this? It leads to better, more reliable, and more reusable code

The `Board` struct

struct Board{T<:Number}
   grid::Matrix{T} 
end
Base.size(b::Board) = size(b.grid)
Base.getindex(b::Board, ci::CartesianIndex) = b.grid[ci]
Base.setindex!(b::Board, val, ci::CartesianIndex) = setindex!(b.grid, val, ci)
Base.eachindex(board::Board) = CartesianIndices(size(board.grid))
Base.similar(b::Board{T}) where T = Board(zeros(T, size(b)))

Base.in(i::CartesianIndex, board::Board) = begin 
    x,y = size(board)
    i[1] <= x && i[1] > 0 && i[2] > 0 && i[2] <= y
end

Board(p = 0.3, sz = (100,100)) = Board(rand(Bernoulli(p), sz))

Board

The simulation loop

function simulate(init::Board, numtimesteps=100)
    timeseries = [similar(init) for _ in 1:numtimesteps]
    timeseries[begin] = init
    
    for t in 2:numtimesteps
        update!(timeseries[t-1], timeseries[t])
    end
    
    return timeseries
end

simulate (generic function with 2 methods)

So we’re done? No, we need to implement the update! function

function update!(oldboard, newboard)
    for cellindex in eachindex(oldboard)
        update!(oldboard, newboard, cellindex)
    end
end

update! (generic function with 1 method)

and now the update! function for each index

 function update!(oldboard, newboard, cellindex)
    st = oldboard[cellindex]
    numneighbors = countneighbors(oldboard, cellindex)
    newboard[cellindex] = nextstate(st, numneighbors)
end

update! (generic function with 2 methods)

Now we need to final functions: countneighbors and nextstate.

function countneighbors(board, index)
    offsets = filter(x->x!=CartesianIndex(0,0), CartesianIndices((-1:1,-1:1)))
    neighborindices = [index + o for o in offsets]
    filter!(x-> x ∈ board, neighborindices)
    return sum(board.grid[neighborindices])
end

countneighbors (generic function with 1 method)

function nextstate(state, numlivingneighbors)
    if !state 
        return numlivingneighbors == 3
    elseif numlivingneighbors == 2 || numlivingneighbors == 3
        return true
    elseif numlivingneighbors < 2 || numlivingneighbors > 3 
        return false
    end
end

nextstate (generic function with 1 method)

@time simulate(Board(0.1,(50,50)), 100);

  0.128976 seconds (990.10 k allocations: 121.109 MiB, 26.04% gc time)

    @time simulate(Board(0.1,(50,50)), 100);

  0.126039 seconds (990.10 k allocations: 121.109 MiB, 28.53% gc time)

function animate(timeseries) 
    pltsettings = (
    aspectratio=1, 
    cbar=:none, 
    frame=:box,
    c=:binary)
    return @animate for t in 1:length(timeseries)
        heatmap(timeseries[t].grid;
            pltsettings..., 
            xlim=(1,size(timeseries[begin])[1]),
            ylim=(1,size(timeseries[begin])[2])
        )
    end
end

animate (generic function with 1 method)

seq = simulate(Board(0.2,(50,50)), 100);
gif(animate(seq), fps=10)

┌ Info: Saved animation to 
│   fn = /home/michael/phd/GOL/tmp.gif
└ @ Plots /home/michael/.julia/packages/Plots/1KWPG/src/animation.jl:114

Adding Patterns

abstract type Pattern end
Base.size(s::Pattern) = s.sz
grid(s::Pattern) = s.shape

grid (generic function with 1 method)

struct Glider <: Pattern 
    sz
    grid
end 
function Glider() 
    mat =  [false false true;
            true  false true;
            false true  true];
    return Glider(size(mat), mat)
end
(g::Glider)() = g.grid

struct Blinker <: Pattern 
    sz
    grid
end 
function Blinker()
    mat = [false true false;
         false true false;
         false true false]
    return Glider(size(mat), mat)
end
(b::Blinker)() = b.grid

struct Pulsar <: Pattern 
    sz
    grid
end 
(p::Pulsar)() = p.grid
function Pulsar() 
    mat =  [0 0 0 0 1 0 0 0 0 0 1 0 0 0 0;
            0 0 0 0 1 0 0 0 0 0 1 0 0 0 0;
            0 0 0 0 1 1 0 0 0 1 1 0 0 0 0;
            0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
            1 1 1 0 0 1 1 0 1 1 0 0 1 1 1;
            0 0 1 0 1 0 1 0 1 0 1 0 1 0 0;
            0 0 0 0 1 1 0 0 0 1 1 0 0 0 0;
            0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
            0 0 0 0 1 1 0 0 0 1 1 0 0 0 0;
            0 0 1 0 1 0 1 0 1 0 1 0 1 0 0;
            1 1 1 0 0 1 1 0 1 1 0 0 1 1 1;
            0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
            0 0 0 0 1 1 0 0 0 1 1 0 0 0 0;
            0 0 0 0 1 0 0 0 0 0 1 0 0 0 0;
            0 0 0 0 1 0 0 0 0 0 1 0 0 0 0]
    return Pulsar(size(mat), mat)
end

Pulsar

# Note: no boundscheck

function add!(board::Board, shape::T, base::Tuple) where T<:Pattern
    sz = size(shape)
    upper = base .+ sz .- 1
    board.grid[base[1]:upper[1],base[2]:upper[2]] .= shape()
    board
end

add! (generic function with 1 method)

add!(Board(), Glider(), (10,10))

Board{Bool}(Bool[0 0 … 0 0; 0 0 … 0 0; … ; 0 1 … 1 0; 0 1 … 0 0])

Base.:+(board::Board, s::Shape) = begin
    add!(board, s)
    return board
end

board = Board(zeros(Bool, 70,70))
add!(board, Glider(), (10,10))
add!(board, Blinker(), (35,50))
add!(board, Pulsar(), (30,30))

Board{Bool}(Bool[0 0 … 0 0; 0 0 … 0 0; … ; 0 0 … 0 0; 0 0 … 0 0])

heatmap(board.grid, c=:binary, aspectratio=1);

out = simulate(board, 500);

gif(animate(out), fps=10)

┌ Info: Saved animation to 
│   fn = /home/michael/phd/GOL/tmp.gif
└ @ Plots /home/michael/.julia/packages/Plots/1KWPG/src/animation.jl:114

DynamicGrids.jl

GOL

init = rand(Bernoulli(0.3), (100,100));

@time out = sim!(ArrayOutput(init, tspan=1:100), Life());

  0.024749 seconds (561 allocations: 1.036 MiB, 72.43% compilation time)

const DEAD, ALIVE, BURNING = 1, 2, 3

neighbors_rule = let prob_combustion=0.0001, prob_regrowth=0.01
    Neighbors(Moore(1)) do data, neighborhood, cell, I
        if cell == ALIVE
            if BURNING in neighborhood
                BURNING
            else
                rand() <= prob_combustion ? BURNING : ALIVE
            end
        elseif cell == BURNING
            DEAD
        else
            rand() <= prob_regrowth ? ALIVE : DEAD
        end
    end
end

[31mNeighbors[39m(
    f = var"#9#10"{Float64, Float64}
    neighborhood = Moore{1, 2, 8, Nothing}
)

init = fill(ALIVE, 200, 200)
output = ArrayOutput(init; 
    tspan=1:500, 
    fps=15, 
);

@time sim!(output, neighbors_rule);

  1.356233 seconds (3.26 M allocations: 177.713 MiB, 11.26% gc time, 73.44% compilation time)

cs = cgrad([:black, :green, :orange]);

anim = @animate for i in 1:length(output)
    heatmap(output[i], c=cs, clim=(1,3))
end

Animation("/tmp/jl_NczJPE", ["000001.png", "000002.png", "000003.png", "000004.png", "000005.png", "000006.png", "000007.png", "000008.png", "000009.png", "000010.png"  …  "000491.png", "000492.png", "000493.png", "000494.png", "000495.png", "000496.png", "000497.png", "000498.png", "000499.png", "000500.png"])

gif(anim)

┌ Info: Saved animation to 
│   fn = /home/michael/phd/GOL/tmp.gif
└ @ Plots /home/michael/.julia/packages/Plots/1KWPG/src/animation.jl:114